ts_data3 <- ts(read.csv("ts3.txt", header = TRUE, as.is = FALSE, sep = ','))
plot(ts_data3)
ts_data7 <- ts(read.csv("ts7.txt", header = TRUE, as.is = FALSE, sep = ','))
plot(ts_data7)
Ряд 3 не колеблется вокруг 0, значит в модели есть константа, ряд нестационарный.
Ряд 7 колеблется вокруг 0, значит в модели нет константы, предположим, что ряд стационарный.
acf(ts_data3)
pacf(ts_data3)
adf.test(ts_data3)
## Augmented Dickey-Fuller Test
## alternative: stationary
##
## Type 1: no drift no trend
## lag ADF p.value
## [1,] 0 0.7251 0.853
## [2,] 1 0.1004 0.673
## [3,] 2 0.0881 0.669
## [4,] 3 0.0717 0.665
## [5,] 4 0.0314 0.653
## [6,] 5 0.0584 0.661
## [7,] 6 0.0642 0.662
## Type 2: with drift no trend
## lag ADF p.value
## [1,] 0 -2.70 0.0788
## [2,] 1 -2.07 0.3017
## [3,] 2 -2.05 0.3096
## [4,] 3 -2.11 0.2826
## [5,] 4 -2.06 0.3034
## [6,] 5 -2.19 0.2530
## [7,] 6 -2.20 0.2492
## Type 3: with drift and trend
## lag ADF p.value
## [1,] 0 -3.20 0.0876
## [2,] 1 -2.29 0.4550
## [3,] 2 -2.27 0.4614
## [4,] 3 -2.34 0.4338
## [5,] 4 -2.28 0.4593
## [6,] 5 -2.42 0.4007
## [7,] 6 -2.44 0.3923
## ----
## Note: in fact, p.value = 0.01 means p.value <= 0.01
Смотрим на p-value для модели с константой (2 и 3), p-value во всех случаях больше 0.05, считаем, что присутствует единичный корень, а значит продифференцируем ряд.
ts_data3_diff <- diff(ts_data3, differences = 1)
plot(ts_data3_diff)
adf.test(ts_data3_diff)
## Augmented Dickey-Fuller Test
## alternative: stationary
##
## Type 1: no drift no trend
## lag ADF p.value
## [1,] 0 -12.73 0.01
## [2,] 1 -11.70 0.01
## [3,] 2 -10.39 0.01
## [4,] 3 -9.51 0.01
## [5,] 4 -8.98 0.01
## [6,] 5 -8.69 0.01
## [7,] 6 -8.07 0.01
## Type 2: with drift no trend
## lag ADF p.value
## [1,] 0 -12.73 0.01
## [2,] 1 -11.70 0.01
## [3,] 2 -10.39 0.01
## [4,] 3 -9.51 0.01
## [5,] 4 -8.99 0.01
## [6,] 5 -8.69 0.01
## [7,] 6 -8.07 0.01
## Type 3: with drift and trend
## lag ADF p.value
## [1,] 0 -12.98 0.01
## [2,] 1 -11.96 0.01
## [3,] 2 -10.65 0.01
## [4,] 3 -9.75 0.01
## [5,] 4 -9.26 0.01
## [6,] 5 -8.98 0.01
## [7,] 6 -8.36 0.01
## ----
## Note: in fact, p.value = 0.01 means p.value <= 0.01
Теперь ряд колеблется вокруг 0. И при этом p-value в модели без сдвига случаях меньше 0.05, считаем, что нет единичного корня и \(d\) = 1, так как хватило одного дифференцирования.
acf(ts_data3_diff)
pacf(ts_data3_diff)
На графике acf видим убывание, похожее на \(e^x\) при этом на графике pacf значима только первая корреляция, считаем \(p\) = 1 и \(q\) = 0. (Если было бы \(q > 0\), то acf быстрее бы убывало, при этом \(p\) определяем как количество значимых корреляций на pacf)
Получили модель ARIMA(1,1,0). Положительное \(\varphi \approx 0.7\).
aa3 <- auto.arima(ts_data3, test = 'adf')
aa3
## Series: ts_data3
## ARIMA(2,1,1)
##
## Coefficients:
## ar1 ar2 ma1
## -0.2237 0.6776 0.9585
## s.e. 0.0713 0.0578 0.0610
##
## sigma^2 = 1.068: log likelihood = -1450.93
## AIC=2909.86 AICc=2909.9 BIC=2929.49
Выходит совсем другая модель :)
model3_my <- Arima(ts_data3, order = c(1, 1, 0))
checkresiduals(model3_my)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,0)
## Q* = 8.2748, df = 9, p-value = 0.5067
##
## Model df: 1. Total lags used: 10
checkresiduals(aa3)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,1,1)
## Q* = 6.574, df = 7, p-value = 0.4745
##
## Model df: 3. Total lags used: 10
В обеих случаях остаток- белый шум \(\sim N(0,\sigma^2)\).
Построим прогноз для обеих моделей (но мне больше нравится моя модель, так как в ней меньше параметров).
Отрежем от ряда 65 точек (после этой точки резкое убывание) и построим прогнозы по моделям, подогнанным под обрезанный ряд.
ts_data3_fc <- ts_data3[1:935]
plot(forecast::forecast(Arima(ts_data3_fc, order = c(1, 1, 0)), h = 65))
lines(936:1000, ts_data3[936:1000], col = 'red')
plot(forecast::forecast(Arima(ts_data3_fc, order = c(2, 1, 1)), h = 65))
lines(936:1000, ts_data3[936:1000], col = 'red')
Прогнозы очень близки и при этом из-за резкого убывания графика рельное продолжение попадает только в 95% интервал.
adf.test(ts_data7)
## Augmented Dickey-Fuller Test
## alternative: stationary
##
## Type 1: no drift no trend
## lag ADF p.value
## [1,] 0 -7.50 0.01
## [2,] 1 -8.76 0.01
## [3,] 2 -8.12 0.01
## [4,] 3 -7.59 0.01
## [5,] 4 -7.72 0.01
## [6,] 5 -7.49 0.01
## Type 2: with drift no trend
## lag ADF p.value
## [1,] 0 -7.58 0.01
## [2,] 1 -8.88 0.01
## [3,] 2 -8.24 0.01
## [4,] 3 -7.72 0.01
## [5,] 4 -7.88 0.01
## [6,] 5 -7.68 0.01
## Type 3: with drift and trend
## lag ADF p.value
## [1,] 0 -7.60 0.01
## [2,] 1 -8.91 0.01
## [3,] 2 -8.27 0.01
## [4,] 3 -7.75 0.01
## [5,] 4 -7.91 0.01
## [6,] 5 -7.70 0.01
## ----
## Note: in fact, p.value = 0.01 means p.value <= 0.01
p-value для модели без константы во всех случаях меньше 0.05, не принимаем гипотезу о наличии единичного корня (процесс стационарный). В модели считаем \(d\) = 0.
acf(ts_data7)
pacf(ts_data7)
На графике acf видим убывание, похожее на \(e^x\cdot\cos{x}\) и на графике pacf значима две корреляции, считаем \(p\) = 2 и \(q\) = 0.
Получили модель ARIMA(2,0,0).
aa7 <- auto.arima(ts_data7, test = 'adf')
aa7
## Series: ts_data7
## ARIMA(2,0,0) with zero mean
##
## Coefficients:
## ar1 ar2
## 0.9625 -0.2080
## s.e. 0.0437 0.0438
##
## sigma^2 = 0.9723: log likelihood = -701.99
## AIC=1409.97 AICc=1410.02 BIC=1422.61
checkresiduals(aa7)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,0,0) with zero mean
## Q* = 3.48, df = 8, p-value = 0.9007
##
## Model df: 2. Total lags used: 10
Вау, совпало.
Отрежем от ряда 31 точку (после этой точки резкое возрастание) и построим прогноз по модели, подогнанной под обрезанный ряд.
ts_data7_fc <- ts_data7[1:469]
plot(forecast::forecast(Arima(ts_data7_fc, order = c(2, 0, 0)), h = 31))
lines(470:500, ts_data7[470:500], col = 'red')
Реальное продолжение ряда помещается в 95% интервал и колеблется вокруг спрогнозированной линии.
ts_data <- read.csv("MRTSSM4453USN.csv", header = TRUE, as.is = FALSE, sep = ',')[,1:2]
ts_data[,1] <- c(0:359)
ts1 <- ts_data[1:336,]
ts1 <- ts(ts1[,2], start = c(1992, 1), frequency = 12)
ts1s <- ts(ts1[1:324], start = c(1992, 1), frequency = 12)
ts1s2 <- ts(ts1[1:312], start = c(1992, 1), frequency = 12)
ssa.result_1 <- ssa(ts1s, L=324/2)
ssa.result_2 <- ssa(ts1s2, L=312/2)
Для одного периода:
auto.arima(ts1s)
## Series: ts1s
## ARIMA(3,1,2)(0,1,2)[12]
##
## Coefficients:
## ar1 ar2 ar3 ma1 ma2 sma1 sma2
## -0.3078 0.1964 0.3714 -0.5782 -0.3112 -0.2548 -0.1585
## s.e. 0.1799 0.0944 0.0719 0.1854 0.1438 0.0608 0.0530
##
## sigma^2 = 5872: log likelihood = -1788.86
## AIC=3593.73 AICc=3594.2 BIC=3623.65
bf_v <- forecast::forecast(ssa.result_1, groups = list(1:15), method = "vector", bootstrap = TRUE, len = 12, R = 100, level = 0.95, interval = 'prediction')
plot(bf_v)
lines(ts1, col = 'black', t='l')
bvec_1 <- rmse(ts1[325:336], bf_v$mean)
ar_1 <- forecast::forecast(Arima(ts1s, order = c(3, 1, 2), seasonal = c(0, 1, 2)), h = 12)
plot(ar_1)
lines(ts1, col = 'black', t='l')
fora_1 <- rmse(ts1[325:336], ar_1$mean)
ets_1 <- forecast::forecast(ets(ts1s), h = 12)
forets_1 <- rmse(ts1[325:336], ets_1$mean)
plot(ets_1)
lines(ts1, col = 'black', t='l')
bvec_1
## [1] 85.11311
fora_1
## [1] 97.87564
forets_1
## [1] 120.1644
Лучший результат– SSA с MSE = 85.11.
auto.arima(ts1s2)
## Series: ts1s2
## ARIMA(3,1,2)(0,1,2)[12]
##
## Coefficients:
## ar1 ar2 ar3 ma1 ma2 sma1 sma2
## -0.2954 0.2044 0.3700 -0.5788 -0.3164 -0.2615 -0.1522
## s.e. 0.1867 0.0972 0.0742 0.1921 0.1480 0.0632 0.0546
##
## sigma^2 = 5792: log likelihood = -1717.72
## AIC=3451.44 AICc=3451.93 BIC=3481.04
bf_r <- forecast::forecast(ssa.result_2, groups = list(1:12), method = "recurrent", bootstrap = TRUE, len = 24, R = 100, level = 0.95, interval = 'prediction')
plot(bf_r)
lines(ts1, col = 'black', t='l')
brec_2 <- rmse(ts1[313:336], bf_r$mean)
ar_2 <- forecast::forecast(Arima(ts1s2, order = c(3, 1, 2), seasonal = c(0, 1, 2)), h = 24)
plot(ar_2)
lines(ts1, col = 'black', t='l')
fora_2 <- rmse(ts1[313:336], ar_2$mean)
ets_2 <- forecast::forecast(ets(ts1s2), h = 24)
forets_2 <- rmse(ts1[313:336], ets_2$mean)
plot(ets_2)
lines(ts1, col = 'black', t='l')
brec_2
## [1] 134.0593
fora_2
## [1] 107.2831
forets_2
## [1] 129.947
Лучший результат– seasonal ARIMA с MSE = 107.28.
Построим два модельных ряда, первый– квадратичный тренд + синус с большим периодом + \(e^x\cdot cos(2\pi\omega n)\) + достаточно сильный шум; второй– квадратичный тренд + \(e^x\cdot cos(2\pi\omega n)\) с меньшей амплитудой + шум слабее.
a <- -1/600
b <- 12
c <- 5
n <- c(1:1008)
ts_mod1 <- ts(c*exp(a*n)*cos(2/b*pi*n) + (n/200)^2 + 100 + sin(2*pi*n/500)) + rnorm(1008, 0, 0.45)
ts_mod2 <- ts(2*c/3*exp(a*n)*cos(2/b*pi*n) + (n/200)^2 + 100) + rnorm(1008, 0, 0.15)
plot(ts_mod1, type = 'l')
plot(ts_mod1, type = 'l')
lines(ts_mod2, type = 'l', col = 'red')
Сравним результаты ssa и mssa:
ssa_signal_1 <- ssa(ts_mod1, svd.method = 'svd')
plot(ssa_signal_1, type = 'vectors')
ssa_signal_2 <- ssa(ts_mod2, svd.method = 'svd')
plot(ssa_signal_2, type = 'vectors')
ssa_signal_m <- ssa(cbind(ts_mod1, ts_mod2), svd.method = 'svd', kind = 'mssa')
plot(ssa_signal_m, type = 'vectors')
r1 <- reconstruct(ssa_signal_1,
groups = list(Trend = c(1, 4:6),
Seasonality = c(2:3)))
r2 <- reconstruct(ssa_signal_m,
groups = list(Trend = c(1, 4:6),
Seasonality = c(2:3)))
plot(r2, add.residuals = FALSE,
plot.method = "xyplot",
slice = list(component = 1),
screens = list('ts_mod1', 'ts_mod2'),
col =
c("blue", "green"),
lty = rep(c(1, 2), each = 2),
scales = list(y = list(draw = FALSE)),
layout = c(1, 2))
plot(r2, plot.method = "xyplot", add.original = FALSE,
add.residuals = FALSE, slice = list(component = 2),
col =
c("blue", "green"),
scales = list(y = list(draw = FALSE)),
layout = c(1, 2))
plot(r1$Trend + r1$Seasonality, t = 'l')
plot(r2$Trend[1:1008] + r2$Seasonality[1:1008], t = 'l')
plot(ts_mod1, t = 'l')
lines(r1$Trend, t = 'l', col = 'blue')
lines(r2$Trend[1:1008], t = 'l', col = 'red')
Восстановленные ряды не отличаются, тренды отличны только в начале при этом очень слабо.
На данном примере результаты очень близки, нельзя сказать, что какой-то из методов показывает результаты лучше.
Взяли изображение с шумом.
img <- readJPEG("12.jpg", native =TRUE)
plot(0:2, 0:2, type='n')
rasterImage(img, 0, 0, 2, 2)
2-D SSA с окном (10, 10):
res_ssa <- ssa(img, kind = "2d-ssa", svd.method = "eigen", L = c(10,10))
plot(res_ssa, type = "vectors", idx = 1:20,
cuts = 255, layout = c(10, 2),
plot.contrib = FALSE)
plot(wcor(res_ssa, groups = 1:30),
scales = list(at = c(10, 20, 30)))
r <- reconstruct(res_ssa, groups = list(Signal = c(1:4, 16)))
plot(r, cuts = 255, layout = c(3, 1))
В результате шум удален, возможно не полностью, попробуем другую длину окна (26,26):
res_ssa <- ssa(img, kind = "2d-ssa", svd.method = "eigen", L = c(26,26))
plot(res_ssa, type = "vectors", idx = 1:20,
cuts = 255, layout = c(10, 2),
plot.contrib = FALSE)
plot(wcor(res_ssa, groups = 1:30),
scales = list(at = c(10, 20, 30)))
r <- reconstruct(res_ssa, groups = list(Signal = c(1:10)))
plot(r, cuts = 255, layout = c(3, 1))
Шум удаляется, контуры чуть более четкие, но понадобилось большее число компонент.